qbit=1;showfig = 0;
close all
Nsub=256;   % no. of subcarriers for rate computation
Ngrid = 1000;
L = 1;

% factor = 1.6;
% x1 = 40e-1*factor;
% y1 = 50e-1*ones(size(factor));
% z1 = 90e-1*factor;
% x2 = -30e-1*factor;
% y2 = 50e-1*ones(size(factor));
% z2 = 40e-1*factor;
% TRX = zeros(6,length(factor));
% TRX(1,:) = x1;
% TRX(2,:) = y1;
% TRX(3,:) = z1;
% TRX(4,:) = x2;
% TRX(5,:) = y2;
% TRX(6,:) = z2;
load('TRX_std.mat','TRX','NTRX')
NTRX = 200;
dist = rms(TRX(1:3,:))*sqrt(3)+rms(TRX(4:6,:))*sqrt(3);

fc = 30e9;
c = 3e8;
BW_des = 1.5e9;
Ptx = 0:5:30;%dBm

cal_rate = 0;
if cal_rate

Rates10 = zeros(1,length(Ptx),NTRX);
Rates3 = zeros(1,length(Ptx),NTRX);
Rates2 = zeros(1,length(Ptx),NTRX);
Rates1 = zeros(1,length(Ptx),NTRX);
Rates0 = zeros(3,length(Ptx),NTRX);
Ratessub = zeros(3,length(Ptx),NTRX);
for n = 1:length(Ptx)
    ptx = Ptx(n);
    parfor m = 1:NTRX
    
   Rates10(:,n,m) = compute_rate_GS(10,Nsub,L,TRX(:,m),fc,BW_des,ptx,1); 
    Rates1(:,n,m) = compute_rate_GS(1,Nsub,L,TRX(:,m),fc,BW_des,ptx,showfig); 
    Rates2(:,n,m) = compute_rate_GS(2,Nsub,L,TRX(:,m),fc,BW_des,ptx,showfig); 
    Rates3(:,n,m) = compute_rate_GS(3,Nsub,L,TRX(:,m),fc,BW_des,ptx,showfig); 
    Rates0(:,n,m) = my_compute_rate(10,Nsub,L,TRX(:,m),fc,BW_des,ptx,showfig); 
    Ratessub(1,n,m) = compute_rate_sub(10,Nsub,L,TRX(:,m),fc,BW_des,ptx,showfig,[2 0]);
    end
    disp(n)
end
Rates2 = mean(Rates2,3);
Rates1 = mean(Rates1,3);
Rates3 = mean(Rates3,3);
Rates10 = mean(Rates10,3);
Rates0 = mean(Rates0,3);
Ratessub = mean(Ratessub,3);
save('data/Rates_sum_rate_vs_power.mat','Rates10',"Rates3","Rates1","Rates2",'Rates0','Ptx','Ratessub')
else
    load('data/Rates_sum_rate_vs_power.mat');
end


plot(Ptx,Rates0(3,:),'--','LineWidth',2) % optimal
hold on
grid on
plot(Ptx,Rates10(1,:),'v-','LineWidth',2) % proposed
plot(Ptx,Rates3(1,:),'o-','LineWidth',2) % proposed b = 3
plot(Ptx,Rates2(1,:),'*-.','LineWidth',2) % proposed b = 2
plot(Ptx,Rates1(1,:),'s-','LineWidth',2) % proposed b = 1
plot(Ptx,Rates0(2,:),'v--','LineWidth',2) %FZ-SPM
plot(Ptx,Rates0(1,:),'k^--','LineWidth',2) % classial 
% plot(Ptx,Rates1(1,:),'d-','LineWidth',2) % classial b=1
plot(Ptx,Ratessub(1,:),'+--','LineWidth',2)
h=plot(dist(round(end/2)),0);
h = h.Parent;
xlabel('Transmit power (dBm)')
ylabel('Average achievable rate (bit/s)')
legend('Optimal beamforming','FZ-GSA','FZ-GSA, Q_b = 3' ...
    ,'FZ-GSA, Q_b = 2','FZ-GSA, Q_b = 1','FZ-SPM','Classical beamforming','VSA-based beamforming' ...
    ,'Location','northwest', 'FontSize', 14)
hold off
set(get(gca, 'Legend'), 'Box', 'off'); 
set(gca,'FontName','Times New Roman')
set(gca, 'FontSize', 16);
set(gcf, 'Position', [100 100 800 650]); % 调整当前图窗的大小
set(0, 'DefaultAxesFontName', 'Times New Roman');

saveas(gcf,'fig/Fig_sum_rate_vs_power.fig')
saveas(gcf,'fig/Fig_sum_rate_vs_power.png')